/**
 * @file test_require_b.cpp
 * @brief This is used to satisty the requirement (b) of PROJECT1
 * @author Linhuo (1066117119@qq.com)
 * @version 1.0
 * @date 2021-04-28
 * 
 * @copyright Copyright (c) 2021  Linhuo
 * 
 * @par 修改日志:
 * <table>
 * <tr><th>Date       <th>Version <th>Author  <th>Description
 * <tr><td>2021-04-29 <td>1.0     <td>Linhuo     <td>内容
 * </table>
 */

#include "../dot_h/TimeIntegrator.h"
#include "../dot_h/ABs.h"
#include "../dot_h/AMs.h"
#include "../dot_h/BDFs.h"
#include "../dot_h/erk4.h"

#include <ctime>

clock_t t[100];
int main(int argc, char *argv[])
{
    /**
     * @brief Initializing for (2)
     */
    Eigen::VectorXd u(6);
    u << 0.994, 0.0, 0.0, 0.0, -2.0015851063790825224, 0.0;
    double T = 17.06521656015796;
    double u1_norm = u.block<3, 1>(0, 0).lpNorm<Eigen::Infinity>();
    int N = 24000;
    double dt = T / N;
    Eigen::MatrixXd test_req_b(6, 4);
    Eigen::VectorXd CPU(4);
    Eigen::VectorXd _Error(4);
    Eigen::MatrixXd display_result(4, 5);
    std::string temp;
    AB test_AB[4], test2_AB[12];
    AM test_AM[4], test2_AM[12];
    BDF test_BDF[4], test2_BDF[12];
    ERK4 test_ERK4, test2_ERK4[3];

    ///#INITIAL CONDITION (2)
    std::cout << "#INITIAL CONDITION (2)" << std::endl;
    /**
 * @brief This is AB solution for (2)
 */
    for (int i = 0; i < 4; i++)
    {
        t[2 * i] = clock();
        temp = "AB" + std::to_string(i + 1);
        test_req_b.col(i) = test_AB[i].Create_as_uwish(temp, u, T, dt);
        t[2 * i + 1] = clock();
        CPU(i) = (double)(t[2 * i + 1] - t[2 * i]) / CLOCKS_PER_SEC;
        _Error(i) = (test_req_b.col(i) - u).block<3, 1>(0, 0).lpNorm<Eigen::Infinity>();
    }

    display_result.col(0) << 1, 2, 3, 4;
    display_result.col(1) << N, N, N, N;
    display_result.col(2) = _Error;
    display_result.col(3) = _Error / u1_norm;
    display_result.col(4) = CPU;

    std::cout << "#name:"
              << "AB" << std::endl;
    std::cout << "      order   steps   Error   Re_Error  CPU_Time" << std::endl;
    std::cout << display_result << std::endl;

    /**
 * @brief This is AM solution for (2)
 */
    for (int i = 0; i < 4; i++)
    {
        t[8 + 2 * i] = clock();
        temp = "AM" + std::to_string(i + 2);
        test_req_b.col(i) = test_AM[i].Create_as_uwish(temp, u, T, dt);
        t[8 + 2 * i + 1] = clock();
        CPU(i) = (double)(t[8 + 2 * i + 1] - t[8 + 2 * i]) / CLOCKS_PER_SEC;
        _Error(i) = (test_req_b.col(i) - u).block<3, 1>(0, 0).lpNorm<Eigen::Infinity>();
    }

    display_result.col(0) << 2, 3, 4, 5;
    display_result.col(1) << N, N, N, N;
    display_result.col(2) = _Error;
    display_result.col(3) = _Error / u1_norm;
    display_result.col(4) = CPU;

    std::cout << "#name:"
              << "AM" << std::endl;
    std::cout << "      order   steps   Error   Re_Error  CPU_Time" << std::endl;
    std::cout << display_result << std::endl;

    /**
 * @brief This is BDF solution for (2)
 */
    for (int i = 0; i < 4; i++)
    {
        t[16 + 2 * i] = clock();
        temp = "BDF" + std::to_string(i + 1);
        test_req_b.col(i) = test_BDF[i].Create_as_uwish(temp, u, T, dt);
        t[16 + 2 * i + 1] = clock();
        CPU(i) = (double)(t[16 + 2 * i + 1] - t[16 + 2 * i]) / CLOCKS_PER_SEC;
        _Error(i) = (test_req_b.col(i) - u).block<3, 1>(0, 0).lpNorm<Eigen::Infinity>();
    }

    display_result.col(0) << 1, 2, 3, 4;
    display_result.col(1) << N, N, N, N;
    display_result.col(2) = _Error;
    display_result.col(3) = _Error / u1_norm;
    display_result.col(4) = CPU;

    std::cout << "#name:"
              << "BDF" << std::endl;
    std::cout << "      order   steps   Error   Re_Error  CPU_Time" << std::endl;
    std::cout << display_result << std::endl;

    /**
 * @brief This is ERK4 solution for (2)
 */
    t[24] = clock();
    temp = "ERK" + std::to_string(4);
    test_req_b.col(0) = test_ERK4.Create_as_uwish(temp, u, T, dt);
    t[25] = clock();
    CPU(0) = (double)(t[25] - t[24]) / CLOCKS_PER_SEC;
    _Error(0) = (test_req_b.col(0) - u).block<3, 1>(0, 0).lpNorm<Eigen::Infinity>();

    display_result.row(0) << 4, N, _Error(0), _Error(0) / u1_norm, CPU(0);
    std::cout << "#name:"
              << "ERK" << std::endl;
    std::cout << "      order   steps   Error   Re_Error  CPU_Time" << std::endl;
    std::cout << display_result.row(0) << std::endl
              << std::endl;

    /**
 * @brief Initializing for (3)
 */
    u << 0.87978, 0, 0, 0, -0.3797, 0;
    T = 19.14045706162071;
    dt = T / N;
    u1_norm = u.block<3, 1>(0, 0).lpNorm<Eigen::Infinity>();
    Eigen::MatrixXd test_temp_2(6, 1);
    Eigen::MatrixXd test_temp_4(6, 1);
    Eigen::MatrixXd test_temp_Explo(6, 3);
    Eigen::VectorXd p_approxi(4);

    std::cout << "#INITIAL CONDITION (3)" << std::endl;
    N = 20000;
    /**
 * @brief This is AB solution for (3).
 */
    for (int i = 0; i < 4; i++)
    {
        if (i == 3)
        {
            dt = T / 1000;
        }
        else if (i == 2)
        {
            dt = T / 6000;
        }
        t[2 * i] = clock();
        temp = "AB" + std::to_string(i + 1);
        test_req_b.col(i) = test2_AB[3 * i].Create_as_uwish(temp, u, T, dt);
        t[2 * i + 1] = clock();
        CPU(i) = (double)(t[2 * i + 1] - t[2 * i]) / CLOCKS_PER_SEC;
        test_temp_2 = test2_AB[3 * i + 1].Create_as_uwish(temp, u, T, dt / 2.0);
        test_temp_4 = test2_AB[3 * i + 2].Create_as_uwish(temp, u, T, dt / 4.0);
        test_temp_Explo.col(0) = (4.0 * test_temp_2 - test_req_b.col(i)) / 3.0;
        test_temp_Explo.col(1) = (4.0 * test_temp_4 - test_temp_2) / 3.0;
        test_temp_Explo.col(2) = (8.0 * test_temp_Explo.col(1) - test_temp_Explo.col(0)) / 7.0;
        _Error(i) = (test_temp_4 - test_temp_Explo.col(2)).block<3, 1>(0, 0).lpNorm<Eigen::Infinity>();
        p_approxi(i) = log2((test_req_b.col(i) - test_temp_2).lpNorm<Eigen::Infinity>()) - log2((test_temp_2 - test_temp_4).lpNorm<Eigen::Infinity>());
    }

    display_result.col(0) << 1, 2, 3, 4;
    display_result.col(1) << N, N, 6000, 1000;
    display_result.col(2) = _Error;
    display_result.col(3) = p_approxi;
    display_result.col(4) = CPU;

    std::cout << "#name:"
              << "AB" << std::endl;
    std::cout << "      order   steps   Error_Approxi   Converg_Rate  CPU_Time" << std::endl;
    std::cout << display_result << std::endl;

    N = 20000;
    /**
 * @brief This is AM solution for (3).
 */
    for (int i = 0; i < 4; i++)
    {
        if (i == 3)
        {
            dt = T / 1000;
        }
        else if (i == 2)
        {
            dt = T / 1000;
        }
        t[2 * i] = clock();
        temp = "AM" + std::to_string(i + 2);
        test_req_b.col(i) = test2_AM[3 * i].Create_as_uwish(temp, u, T, dt);
        t[2 * i + 1] = clock();
        CPU(i) = (double)(t[2 * i + 1] - t[2 * i]) / CLOCKS_PER_SEC;
        test_temp_2 = test2_AM[3 * i + 1].Create_as_uwish(temp, u, T, dt / 2.0);
        test_temp_4 = test2_AM[3 * i + 2].Create_as_uwish(temp, u, T, dt / 4.0);
        test_temp_Explo.col(0) = (4.0 * test_temp_2 - test_req_b.col(i)) / 3.0;
        test_temp_Explo.col(1) = (4.0 * test_temp_4 - test_temp_2) / 3.0;
        test_temp_Explo.col(2) = (8.0 * test_temp_Explo.col(1) - test_temp_Explo.col(0)) / 7.0;
        _Error(i) = (test_temp_4 - test_temp_Explo.col(2)).block<3, 1>(0, 0).lpNorm<Eigen::Infinity>();
        p_approxi(i) = log2((test_req_b.col(i) - test_temp_2).lpNorm<Eigen::Infinity>()) - log2((test_temp_2 - test_temp_4).lpNorm<Eigen::Infinity>());
    }

    display_result.col(0) << 2, 3, 4, 5;
    display_result.col(1) << N, N, 1000, 1000;
    display_result.col(2) = _Error;
    display_result.col(3) = p_approxi;
    display_result.col(4) = CPU;

    std::cout << "#name:"
              << "AM" << std::endl;
    std::cout << "      order   steps   Error_Approxi   Converg_Rate  CPU_Time" << std::endl;
    std::cout << display_result << std::endl;

    N = 30000;
    dt = T / N;
    std::cout << "//loading...(Need half a minute to compute)" << std::endl;
    /**
 * @brief This is BDF solution for (3).
 */
    for (int i = 0; i < 4; i++)
    {
        if (i == 3)
        {
            dt = T / 600;
        }
        else if (i == 2)
        {
            dt = T / 600;
        }
        else if (i == 1)
        {
            dt = T / 2000;
        }
        t[2 * i] = clock();
        temp = "BDF" + std::to_string(i + 1);
        test_req_b.col(i) = test2_BDF[3 * i].Create_as_uwish(temp, u, T, dt);
        t[2 * i + 1] = clock();
        CPU(i) = (double)(t[2 * i + 1] - t[2 * i]) / CLOCKS_PER_SEC;
        test_temp_2 = test2_BDF[3 * i + 1].Create_as_uwish(temp, u, T, dt / 2.0);
        test_temp_4 = test2_BDF[3 * i + 2].Create_as_uwish(temp, u, T, dt / 4.0);
        test_temp_Explo.col(0) = (4.0 * test_temp_2 - test_req_b.col(i)) / 3.0;
        test_temp_Explo.col(1) = (4.0 * test_temp_4 - test_temp_2) / 3.0;
        test_temp_Explo.col(2) = (8.0 * test_temp_Explo.col(1) - test_temp_Explo.col(0)) / 7.0;
        _Error(i) = (test_temp_4 - test_temp_Explo.col(2)).block<3, 1>(0, 0).lpNorm<Eigen::Infinity>();
        p_approxi(i) = log2((test_req_b.col(i) - test_temp_2).lpNorm<Eigen::Infinity>()) - log2((test_temp_2 - test_temp_4).lpNorm<Eigen::Infinity>());
    }

    display_result.col(0) << 1, 2, 3, 4;
    display_result.col(1) << N, 2000, 600, 600;
    display_result.col(2) = _Error;
    display_result.col(3) = p_approxi;
    display_result.col(4) = CPU;

    std::cout << "#name:"
              << "BDF" << std::endl;
    std::cout << "      order   steps   Error_Approxi   Converg_Rate  CPU_Time" << std::endl;
    std::cout << display_result << std::endl;

    N = 1000;
    dt = T / N;
    /**
 * @brief This is ERK4 solution for (3)
 */
    temp = "ERK" + std::to_string(4);
    t[24] = clock();
    test_req_b.col(0) = test2_ERK4[0].Create_as_uwish(temp, u, T, dt);
    t[25] = clock();
    CPU(0) = (double)(t[25] - t[24]) / CLOCKS_PER_SEC;
    test_temp_2 = test2_ERK4[1].Create_as_uwish(temp, u, T, dt / 2.0);
    test_temp_4 = test2_ERK4[2].Create_as_uwish(temp, u, T, dt / 4.0);
    test_temp_Explo.col(0) = (4.0 * test_temp_2 - test_req_b.col(0)) / 3.0;
    test_temp_Explo.col(1) = (4.0 * test_temp_4 - test_temp_2) / 3.0;
    test_temp_Explo.col(2) = (8.0 * test_temp_Explo.col(1) - test_temp_Explo.col(0)) / 7.0;
    _Error(0) = (test_temp_4 - test_temp_Explo.col(2)).block<3, 1>(0, 0).lpNorm<Eigen::Infinity>();
    p_approxi(0) = log2((test_req_b.col(0) - test_temp_2).lpNorm<Eigen::Infinity>()) - log2((test_temp_2 - test_temp_4).lpNorm<Eigen::Infinity>());

    display_result.row(0) << 4, N, _Error(0), p_approxi(0), CPU(0);

    std::cout << "#name:"
              << "ERK" << std::endl;
    std::cout << "      order   steps   Error   Converg_Rate  CPU_Time" << std::endl;
    std::cout << display_result.row(0) << std::endl
              << std::endl;

    std::cout << "Complete!" << std::endl;

    return 0;
};